%% Table 5
% based on estimation results in columns 1 and 2 of Table 4
clc; 
clear;
%-----------------------Parameters--------------------------------------------------------
beta=0.97;
sigma2=0.00172; % variance excess return
numb=100000; % number of simulations
k=1; % set equal to 1 if OLS (column 1 of Table 4), equal to 2 if IV (column 2 of Table 4)
%----------------------END PARAMETERS-----------------------------------------------------
mat=zeros(1,3); % estimated coefficients 
if k==1
  mat=[0.93 0.423 8.243]; % OLS point estimates, Table 4, column 1
  mat1=zeros(3,3); % covariance matrix coefficients Table 4, column 1
  mat1=[0.00001816	0.00002068	-0.00117251
        0.00002068	 0.00074514	-0.00542979
        -0.00117251	-0.00542979	0.66726965]; % OLS
end   

if k==2
  mat=[0.948 0.293 9.403]; % IV point estimates, Table 4, column 2
  mat1=zeros(3,3); % covariance matrix coefficients Table 4, column 2
  mat1=[0.00015956	0.00032383	-0.02909301
        0.00032383	 0.00429428	-0.07169163
        -0.02909301	-0.07169163	7.517093]; % IV
end

a1=mat(1,1); % coefficient z_{i,n,t-1}
a2=mat(1,2); % coefficient val_{i,n,t}
a3=mat(1,3); % coefficient ER_{i,n,t}
delta=a1*beta;

lambda1=(a1-a2)/(a3*(1-delta));
lambda2=a2/(a3*(1-delta));
gamma=(1-a1)/(a3*sigma2);
diff=lambda1-lambda2; 

% next draw "numb" times from the distribution of the estimated parameters
% recompute lambda1, lambda2, gamma and lambda1-lambda2 (diff) each time
eps=normrnd(0,1,3,numb); % 3 by "numb" matrix of random N(0,1) numbers
lab1=zeros(numb,1); %lambda1
lab2=zeros(numb,1); %lambda2
gam=zeros(numb,1); %gamma
dif=zeros(numb,1); %lambda1-lambda2
m=mat1^0.5;
for i=1:numb
    eps1=m*eps(1:3,i);
    a1=mat(1,1)+eps1(1); % new coefficient a1
    a2=mat(1,2)+eps1(2); % new coefficient a2
    a3=mat(1,3)+eps1(3); % new coefficient a3
    lab1(i)=(a1-a2)/(a3*(1-delta));
    lab2(i)=a2/(a3*(1-delta));
    gam(i)=(1-a1)/(a3*sigma2);
    dif(i)=lab1(i)-lab2(i);
end
% sort the values of lambda1, lambda2, gamma and lambda1-lambda2
slambda1=sort(lab1(1:numb)); 
slambda2=sort(lab2(1:numb)); 
sgamma=sort(gam(1:numb)); 
sdiff=sort(dif(1:numb)); 
% compute 2.5 and 97.5 percentiles
lambda1low=slambda1(numb*0.025);
lambda1high=slambda1(numb*0.975);
lambda2low=slambda2(numb*0.025);
lambda2high=slambda2(numb*0.975);
gammalow=sgamma(numb*0.025);
gammahigh=sgamma(numb*0.975);
difflow=sdiff(numb*0.025);
diffhigh=sdiff(numb*0.975);

fprintf('results: point estimate and confidence interval \n\n');
fprintf('lambda1: %6.4f %6.4f %6.4f \n',lambda1,lambda1low,lambda1high)
fprintf('lambda2: %6.4f %6.4f %6.4f \n',lambda2,lambda2low,lambda2high)
fprintf('gamma: %6.4f %6.4f %6.4f \n',gamma,gammalow,gammahigh)
fprintf('lambda1-lambda2: %6.4f %6.4f %6.4f \n',diff,difflow,diffhigh)



